home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Tech Arsenal 1
/
Tech Arsenal (Arsenal Computer).ISO
/
tek-01
/
ddj9304.zip
/
STEREO.ZIP
/
LISTING2.C
< prev
Wrap
INI File
|
1993-01-19
|
21KB
|
521 lines
[LISTING2.C]
/* Created by Victor J. Duvanenko 11-05-90
Program to display the Scanning Tunneling Microscope (STM)
collected image data in 'true' 3-D.
-------------------------------------------------------------------*/
#include <math.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#define BOOLEAN int
#define TRUE 1
#define FALSE 0
#define OK 0
#define NOT_OK -1
#define ACCEPT TRUE
#define REJECT FALSE
#define WHITE 1
#define PI 3.1415926
#define INTEROCULAR_DISTANCE 2.5 /* inches */
#define PIXEL_WIDTH 0.021 /* inches */
#define LEFT_EYE 0
#define RIGHT_EYE 1
#define NUM_DIMENSIONS 3
#define NUM_SHADES 256
#define SCREEN_HEIGHT 512
#define SCREEN_WIDTH 640
#define WINDOW_BORDER_WIDTH 3
#define PARALLEL 0
#define PERSPECTIVE 1
#define NUM_LINES 348 /* for an IBM PC Hercules display */
#define MAX_IMAGE_SIZE 200 /* 400 x 400 samples */
#define MALLOC(x) ((x *) malloc( sizeof(x)))
#define ROUND(x) ((x) > 0.0 ? (int)((x) + 0.5) : (int)((x) - 0.5))
#define MAX(A, B) ((A) > (B) ? (A) : (B))
#define MIN(A, B) ((A) < (B) ? (A) : (B))
#define DEBUG FALSE /* en(dis)able debug sections */
#define DRAW TRUE
#define LEX TRUE
#define PC FALSE
struct point_3D
{
float x, y, z;
};
typedef struct point_3D point_3D_t;
struct point_3D_ex
{
float x, y, z;
float sha, shb; /* shading values for triangle A and B */
};
typedef struct point_3D_ex point_3D_ex_t;
/* STM data file header parameters. */
int num_samples; /* number of samples per scan line */
int log_input; /* 0 - linear input; 1 - logarithmic input */
float z_scale; /* vertical scale */
float scan_sz; /* 0 < X, Y < scan_sz in nm */
point_3D_t min, max; /* min and max of the input/STM data */
/* X, Y, Z values at each sample point - STM gives the Z; X & Y are deduced
from the stepping distance (sampling distance => (scan_sz / num_samples)).*/
point_3D_ex_t image[ MAX_IMAGE_SIZE ][ MAX_IMAGE_SIZE ];
/* Transformed image */
point_3D_ex_t tr_image[ MAX_IMAGE_SIZE ][ MAX_IMAGE_SIZE ];
#if Z_BUFFER
/* Image buffer */
unsigned char b_image[ SCREEN_WIDTH ][ SCREEN_HEIGHT ];
float z_buffer[ SCREEN_WIDTH ][ SCREEN_HEIGHT ];
#endif
/* World window/volume rectangle boundaries - accessable to all routines */
double y_bottom, y_top, x_right, x_left, z_front, z_back;
double y_bottom_d, y_top_d, x_right_d, x_left_d, z_front_d, z_back_d;
point_3D_t light_source; /* light source direction unit vector */
point_3D_t eye_pt; /* viewer eye point in world coordinates */
point_3D_t vpn; /* view plane normal vector */
point_3D_t vup; /* view up vector */
double Tm[ 4 ][ 4 ]; /* transformation matrix */
double proj_plane; /* perspective projection plane on -Z-axis */
/* some global variables */
int f_color;
/* External functions */
extern void set_clip_volume( double, double, double, double, double, double );
extern find_object_center( double *, double *, double *, int );
extern void read_stm_database();
extern void draw_clip_window();
extern matrix_mult( double *, int, int, double *, int, int, double * );
/*-------------------------------------------------------------------
Procedure that creates a composite viewing transformation matrix
in accordance with p. 258 - 261 of 1st ed. of F & vD.
The viewer eye point will be moved to the origin looking down the
negative Z-axis with the view up vector in the yz-plane in positive
y half of the plane (pointing up).
P1 = VRP = eye_pt
P2 = origin
P3 = P1 + vup
This procedure is executed only once and only when the viewer eye
point position changes.
---------------------------------------------------------------------*/
make_composite_viewing_matrix( proj )
int proj; /* PARALLEL or PERSPECTIVE */
{
double p1[4], p2[4], p3[4], p_tmp[4];
double T[4][4], Rx[4][4], Ry[4][4], Rz[4][4], C_tmp1[4][4], C_tmp2[4][4];
double Per[4][4], d1, d2, d12, cos_ang, sin_ang;
/* Initialize the three points */
p1[0] = eye_pt.x; p1[1] = eye_pt.y; p1[2] = eye_pt.z; p1[3] = 1.0;
p2[0] = p2[1] = p2[2] = 0.0; p2[3] = 1.0;
p3[0] = p1[0] + vup.x; p3[1] = p1[1] + vup.y; p3[2] = p1[2] + vup.z;
p3[3] = 1.0;
/* Magnitude of vector p1->p2 */
d12 = sqrt( p1[0] * p1[0] + p1[1] * p1[1] + p1[2] * p1[2] );
/* Create the translation matrix. */
set_to_identity( T, 4 );
T[3][0] = -p1[0]; T[3][1] = -p1[1]; T[3][2] = -p1[2];
/* Translate the three points p1, p2, and p3 to the origin. */
matrix_mult( p1, 1, 4, T, 4, 4, p_tmp );
copy_matrix( p_tmp, p1, 1, 4 );
matrix_mult( p2, 1, 4, T, 4, 4, p_tmp );
copy_matrix( p_tmp, p2, 1, 4 );
matrix_mult( p3, 1, 4, T, 4, 4, p_tmp );
copy_matrix( p_tmp, p3, 1, 4 );
d1 = sqrt( p2[0] * p2[0] + p2[2] * p2[2] ); /* length of projection */
cos_ang = -p2[2] / d1;
sin_ang = p2[0] / d1;
/* Create the rotation about Y-axis matrix. */
set_to_identity( Ry, 4 );
Ry[0][0] = cos_ang; Ry[0][2] = -sin_ang;
Ry[2][0] = sin_ang; Ry[2][2] = cos_ang;
/* Rotate the three points p2, and p3 about the Y-axis. */
/* p1 is at the origin after translation => no need to rotate. */
matrix_mult( p2, 1, 4, Ry, 4, 4, p_tmp );
copy_matrix( p_tmp, p2, 1, 4 );
matrix_mult( p3, 1, 4, Ry, 4, 4, p_tmp );
copy_matrix( p_tmp, p3, 1, 4 );
cos_ang = -p2[2] / d12;
sin_ang = -p2[1] / d12;
/* Create the rotation about X-axis matrix. */
set_to_identity( Rx, 4 );
Rx[1][1] = cos_ang; Rx[1][2] = sin_ang;
Rx[2][1] = -sin_ang; Rx[2][2] = cos_ang;
/* Rotate the three points p2, and p3 about the X-axis. */
matrix_mult( p2, 1, 4, Rx, 4, 4, p_tmp );
copy_matrix( p_tmp, p2, 1, 4 );
matrix_mult( p3, 1, 4, Rx, 4, 4, p_tmp );
copy_matrix( p_tmp, p3, 1, 4 );
/* Sanity check. */
printf( "The view vector should be [ %lf %lf %lf %lf ]\n", 0.0, 0.0, -d12,
1.0 );
printf( "The view vector is [ %lf %lf %lf %lf ]\n", p2[0], p2[1],
p2[2], p2[3] );
d2 = sqrt( p3[0] * p3[0] + p3[1] * p3[1] );
cos_ang = p3[1] / d2;
sin_ang = p3[0] / d2;
/* Create the rotation about Z-axis matrix. */
set_to_identity( Rz, 4 );
Rz[0][0] = cos_ang; Rz[0][1] = sin_ang;
Rz[1][0] = -sin_ang; Rz[1][1] = cos_ang;
/* At this point the translation, and all rotation matrices are known
and need to be combined into a single transformaation matrix. */
matrix_mult( T, 4, 4, Ry, 4, 4, C_tmp1 );
matrix_mult( C_tmp1, 4, 4, Rx, 4, 4, C_tmp2 );
matrix_mult( C_tmp2, 4, 4, Rz, 4, 4, C_tmp1 );
copy_matrix( C_tmp1, Tm, 4, 4 );
}
/*-------------------------------------------------------------------
Procedure that prints the main menu of operations available to a user.
---------------------------------------------------------------------*/
void print_menu()
{
printf( "Load STM database ...........................(1)\n" );
printf( "Set light source direction vector ...........(2)\n" );
printf( "Set viewer reference point ..................(3)\n" );
printf( "Set view up vector ..........................(4)\n" );
printf( "Display with parallel projection .........(5)\n" );
printf( "Display with perspective projection .........(6)\n" );
printf( "Set perspective projection plane ............(7)\n" );
printf( "Set window limits ...........................(8)\n" );
printf( "Set viewport limits .........................(9)\n" );
printf( "Display in stereo ..........................(10)\n" );
printf( "Quit........................................(26)\n" );
printf( "? " );
}
/* Procedure to transform the image by a transformation matrix.
Assumes the matrix (m) is a 4x4 matrix. The result is homogenized
(divided by W) as well. The source and destination can be the same
array. */
transform_and_homogenize_image( s, d, tm )
point_3D_ex_t s[][ MAX_IMAGE_SIZE ], d[][ MAX_IMAGE_SIZE ];
double *tm; /* transformation matrix */
{
register i, j;
int num_lines;
double p[4]; /* the point to be transformed */
double t[4], inv_W;
num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
for( i = 0; i < num_lines; i++ )
for( j = 0; j < num_lines; j++ )
{
p[0] = s[i][j].x; p[1] = s[i][j].y;
p[2] = s[i][j].z; p[3] = 1.0;
matrix_mult( p, 1, 4, tm, 4, 4, t );
if ( t[3] != 1.0 ) /* divide by W (homogenize) */
{
inv_W = 1.0 / t[3];
t[0] *= inv_W; t[1] *= inv_W; t[2] *= inv_W;
}
d[i][j].x = t[0]; d[i][j].y = t[1]; d[i][j].z = t[2];
d[i][j].sha = s[i][j].sha; d[i][j].shb = s[i][j].shb;
}
}
/* Procedure to render the transformed image. */
render_image( y )
int y;
{
register i, j;
int num_lines;
short v[6], intensity;
num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
num_lines--;
for( i = 0; i < num_lines; i++ )
for( j = 0; j < num_lines; j++ )
{
v[0] = ROUND( tr_image[ i ][ j ].x );
v[1] = ROUND( tr_image[ i ][ j ].y );
v[2] = ROUND( tr_image[ i + 1 ][ j ].x );
v[3] = ROUND( tr_image[ i + 1 ][ j ].y );
v[4] = ROUND( tr_image[ i + 1 ][ j + 1 ].x );
v[5] = ROUND( tr_image[ i + 1 ][ j + 1 ].y );
/* Render triangle A */
intensity = ROUND( tr_image[ i ][ j ].sha * (double)(NUM_SHADES - 1));
if ( intensity > ( NUM_SHADES - 1 ))
intensity = NUM_SHADES - 1; /* saturate */
if ( intensity < 0 )
{
#if 0
printf( "Triangle A, intensity = %d\n", intensity );
printf( "v11.x = %f, v11.y = %f, v11.z = %f\n",
image[i][j].x, image[i][j].y, image[i][j].z );
printf( "v21.x = %f, v21.y = %f, v21.z = %f\n",
image[i+1][j].x, image[i+1][j].y, image[i+1][j].z );
printf( "v22.x = %f, v22.y = %f, v22.z = %f\n",
image[i+1][j+1].x, image[i+1][j+1].y, image[i+1][j+1].z );
#endif
intensity = 0;
}
if ( clip_to_viewport( v, 6, y ) == ACCEPT )
dspoly( &intensity, v, &6 );
v[2] = ROUND( tr_image[ i ][ j + 1 ].x );
v[3] = ROUND( tr_image[ i ][ j + 1 ].y );
/* Render triangle B */
intensity = ROUND( tr_image[ i ][ j ].shb * (double)( NUM_SHADES - 1));
if ( intensity > ( NUM_SHADES - 1 ))
intensity = NUM_SHADES - 1; /* saturate */
if ( intensity < 0 ) intensity = 0;
if ( clip_to_viewport( v, 6, y ) == ACCEPT )
dspoly( &intensity, v, &6 );
}
}
/* Procedure to display the database with a projection. */
display_3D_data( proj )
int proj;
{
double Tw[4][4], S[4][4], Td[4][4], Ry[4][4], Per[4][4], tmp[4][4];
double Left[4][4], Right[4][4];
double e_v, /* interocular distance mapped back to the view coord. */
e_w; /* interocular distance mapped back to the world coord. */
printf( "Computing normals for shading. " );
compute_normals(); /* and the dot products for each triangle */
printf( "Done.\n" );
/* Perspective projection must use three steps:
1) Compute normals and project,
2) Divide by W (homogenize)
3) Transform into the device coordinates. */
if ( proj == PERSPECTIVE )
{
/* map the physical interocular distance into the view port coordinates. */
e_v = INTEROCULAR_DISTANCE / PIXEL_WIDTH; /* e_v == pixels */
/* map from the viewport coordinate system to the world. */
e_w = e_v / (( x_right_d - x_left_d ) / ( x_right - x_left ));
set_to_identity( Left, 4 );
set_to_identity( Right, 4 );
/* Use the Translate, Project, Translate back model. */
/* Create the Left eye transformation matrix. */
set_to_identity( Tw, 4 ); /* translate the left eye to the origin */
Tw[3][0] = -e_w / 2.0;
matrix_mult( Left, 4, 4, Tw, 4, 4, tmp );
/* Create the perspective projection matrix. */
set_to_identity( Per, 4 );
Per[2][3] = 1.0 / proj_plane; /* 1/d */
Per[3][3] = 0.0;
matrix_mult( tmp, 4, 4, Per, 4, 4, Left );
Tw[3][0] = e_w / 2.0; /* translate back */
matrix_mult( Left, 4, 4, Tw, 4, 4, tmp );
copy_matrix( tmp, Left, 4, 4 );
/* Create the Right eye transformation matrix. */
set_to_identity( Tw, 4 ); /* translate the right eye to the origin */
Tw[3][0] = e_w / 2.0;
matrix_mult( Right, 4, 4, Tw, 4, 4, tmp );
/* Create the perspective projection matrix. */
set_to_identity( Per, 4 );
Per[2][3] = 1.0 / proj_plane; /* 1/d */
Per[3][3] = 0.0;
matrix_mult( tmp, 4, 4, Per, 4, 4, Right );
Tw[3][0] = -e_w / 2.0; /* translate back */
matrix_mult( Right, 4, 4, Tw, 4, 4, tmp );
copy_matrix( tmp, Right, 4, 4 );
#if 0
printf( "Transforming, projecting and homogenizing the image. " );
transform_and_homogenize_image( image, tr_image, Tm );
printf( "Done.\n" );
#endif
}
/* Create the world to device transformation matrix. */
/* Create the translation matrix. Translate only in X and Y. */
set_to_identity( Tw, 4 );
Tw[3][0] = -x_left; Tw[3][1] = -y_bottom; Tw[3][2] = 0.0;
/* Create a uniform scale matrix. */
set_to_identity( S, 4 );
S[0][0] = ( x_right_d - x_left_d ) / ( x_right - x_left );
S[1][1] = ( y_top_d - y_bottom_d ) / ( y_top - y_bottom );
S[2][2] = ( z_back_d - z_front_d ) / ( z_back - z_front );
matrix_mult( Tw, 4, 4, S, 4, 4, tmp );
copy_matrix( tmp, Tw, 4, 4 );
/* Create the translation matrix. Translate only in X and Y. */
set_to_identity( Td, 4 );
Td[3][0] = x_left_d; Td[3][1] = y_bottom_d; Td[3][2] = 0.0;
matrix_mult( Tw, 4, 4, Td, 4, 4, tmp );
copy_matrix( tmp, Tw, 4, 4 );
/* Since the device/screen origin on the LEX/90 is in the upper left
we nee to reflect the Y and translate by the
screen height to place the device origin in the bottom left. */
set_to_identity( Ry, 4 );
Ry[1][1] = -1.0;
matrix_mult( Tw, 4, 4, Ry, 4, 4, tmp );
copy_matrix( tmp, Tw, 4, 4 );
set_to_identity( Td, 4 );
Td[3][1] = SCREEN_HEIGHT;
matrix_mult( Tw, 4, 4, Td, 4, 4, tmp );
copy_matrix( tmp, Tw, 4, 4 );
/* Now, Tw has the world to device/screen transformation matrix. */
if ( proj == PARALLEL )
{
/* Beautiful!!! Perform a single transformation of the image (for
parallel projection). */
printf( "Transforming, projecting and mapping the image onto screen. " );
matrix_mult( Tm, 4, 4, Tw, 4, 4, tmp );
copy_matrix( tmp, Tm, 4, 4 );
show_matrix( Tm, 4, 4 );
transform_and_homogenize_image( image, tr_image, Tm );
printf( "Done.\n" );
printf( "Rendering the image. " );
render_image( LEFT_EYE ); printf( "Done.\n" );
}
if ( proj == PERSPECTIVE )
{
printf( "Transforming, projecting, homogenizing and mapping the " );
printf( "Left image onto screen. " );
matrix_mult( Tm, 4, 4, Left, 4, 4, tmp );
matrix_mult( tmp, 4, 4, Tw, 4, 4, Left );
printf( "Left eye transformation matrix.\n" );
show_matrix( Left, 4, 4 );
transform_and_homogenize_image( image, tr_image, Left );
printf( "Done.\n" );
printf( "Rendering the Left eye image. " );
render_image( LEFT_EYE ); printf( "Done.\n" );
printf( "Transforming, projecting, homogenizing and mapping the " );
printf( "Right image onto screen. " );
matrix_mult( Tm, 4, 4, Right, 4, 4, tmp );
matrix_mult( tmp, 4, 4, Tw, 4, 4, Right );
/* Move the right eye view into the lower half of the buffer. */
set_to_identity( Tw, 4 );
Tw[3][1] = SCREEN_HEIGHT; /* move in device coord */
matrix_mult( Right, 4, 4, Tw, 4, 4, tmp );
copy_matrix( tmp, Right, 4, 4 );
printf( "Right eye transformation matrix.\n" );
show_matrix( Right, 4, 4 );
transform_and_homogenize_image( image, tr_image, Right );
printf( "Done.\n" );
dszom( &0, &512, &1 ); /* look at the lower half */
printf( "Rendering the Right eye image. " );
render_image( RIGHT_EYE ); printf( "Done.\n" );
}
}
/*------------------------- MAIN ------------------------------------*/
main( argc, argv )
int argc;
char **argv;
{
register i;
long int j;
int response, tmp;
BOOLEAN done;
double x_min, x_max, y_min, y_max, z_min, z_max;
double x, y, z;
char str[128]; /* a temporary string buffer */
/* setup some defaults */
f_color = 255;
done = FALSE;
light_source.x = 0.0; light_source.y = 0.0; light_source.z = 1.0;
eye_pt.x = 0.0; eye_pt.y = 0.0; eye_pt.z = 260.0;
vpn.x = -(eye_pt.x); vpn.y = -(eye_pt.y); vpn.z = -(eye_pt.z);
vup.x = 0.0; vup.y = 1.0; vup.z = 0.0;
proj_plane = -1.0;
/* A welcome message */
printf( "\n\n\n STM Data Viewer.\n\n\n" );
printf( " Created by Victor J. Duvanenko\n" );
printf( " 11/05/90\n\n\n\n" );
set_clip_volume( -250.0, 250.0, -250.0, 250.0, 0.0, -500.0 );
set_viewport( 50.0, 562.0, 0.0, 512.0, 0.0, 512.0 );
open_lex();
configure_lex();
dsbuff( &0 ); /* display and manipulate buffer 0 */
while( !done )
{
fflush( stdin );
print_menu();
scanf( "%d", &response );
switch( response)
{
case 1: read_stm_database(); dsclr( &(-1) );
set_clip_volume((double)min.x, (double)max.x, (double)min.y,
(double)max.y, (double)min.z, (double)max.z );
display_stm_database(); break;
case 2: set_light_source(); break;
case 3: set_viewer_eye_pt(); break;
case 4: set_view_up_vector(); break;
case 5: dsclr( &(-1));
make_composite_viewing_matrix( PARALLEL );
display_3D_data( PARALLEL ); break;
case 6: dsclr( &(-1));
make_composite_viewing_matrix( PERSPECTIVE );
display_3D_data( PERSPECTIVE ); break;
case 7: set_per_projection_plane(); break;
case 8: printf( "Enter window limits.\n" );
printf( "x_min x_max y_min y_max z_min z_max\n " );
scanf( "%lf%lf%lf%lf%lf%lf", &x_min, &x_max,
&y_min, &y_max, &z_min, &z_max );
set_clip_volume( x_min, x_max, y_min, y_max, z_min,z_max);
break;
case 9: printf( "Enter viewport limits.\n" );
printf( "x_min x_max y_min y_max z_min z_max\n " );
scanf( "%lf%lf%lf%lf%lf%lf", &x_min, &x_max,
&y_min, &y_max, &z_min, &z_max );
set_viewport( x_min, x_max, y_min, y_max,
z_min, z_max );
break;
case 10: dsmov( &512, &640, &2, &1 );
printf( "Type any number to stop the stereo display: " );
scanf( "%d", &tmp );
dszom( &0, &0, &1 ); /* back to one eye view */
break;
case 26: done = TRUE; break;
default : printf( "Illegal choice.\n" );
}
if ( done ) break; /* leave without timing */
}
dscls(); /* close the communication channel to the LEX */
return(0);
}